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HYBRID-DUAL-FOURIER T OMOGRAPHIC 
ALGORITHM FOR A FAST 

THREE-DIMENSIONIAL OPTICAL IMAGE 
RECONSTRUCTION IN TURBID MEDL4 

5 
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This invention was made, in part, with Government 
support awarded by National Aeronautics and Space Admin- 
istration (NASA), and the US Army Medical Research and 
Materiel Command (USAMRMC). The Government may 
have certain rights in this invention. 20 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 

The present invention teaches a novel hybrid-dual -Fourier 25 
inverse algorithm for fast three dimensional (3D) tomo- 
graphic image reconstruction of objects in highly scattering 
media using measured data from multiple sources and 
multiple detectors. This algorithm can be used for noninva- 
sive screening, detection, and diagnosis of cancerous breast 30 
and prostate lesions, and for locating hidden objects in high 
scattering media, such as planes, tanks, objects within an 
animal or human body, or through cloud, fog, or smoke, 
mines under turbid water, and corrosion under paint. 

2. Description of the Related Art 35 

Turbid media blur and make objects inside difficult to be 

detected. Light propagation is diffusive, making the object 
inside not detectable using conventional imaging methods. 

In highly scattering media, when the object is located inside 
a turbid medium with depth greater than 10 scattering 40 
length, the object cannot be seen easily using ballistic light. 

It is also difficult to use time-gated transillumination tech- 
nique to image objects in turbid media with L/I^>20, with L 
the size of the medium and l s the scattering length, because 
of photon starvation of the ballistic light. To overcome this 45 
difficulty one uses the image reconstruction of light in the 
medium to get three-dimensional image. This requires 
understanding of how light travels in the turbid media, and 
an appropriate inverse algorithm. Objects, such as tumors, 
aircraft, and corrosion in highly scattering media can be 50 
imaged using the novel image algorithm in turbid media. 
Turbid media include human tissue, cloud, and under paint. 

Early detection and diagnosis of breast and prostate 
cancers is essential for effective treatment. X-ray mammog- 
raphy, the modality commonly used for breast cancer screen- 55 
ing, cannot distinguish between malignant and benign 
tumors, and is less effective for younger women with dense 
fibrous breasts. If a tumor is suspected from a x-ray mam- 
mogram, a biopsy that requires invasive removal of tissue 
from the suspect region need be performed to determine if 60 
the tumor is benign or malignant. In a majority of the cases, 
the biopsy turns out to be negative, meaning the tumor is 
benign. Besides being subject to an invasive procedure, one 
has to wait an agonizing period until the biopsy results are 
known. A breast cancer screening modality that does not 65 
require tissue removal, and can provide diagnostic informa- 
tion is much desired. 


2 

Prostate cancer has a high incidence of mortality for men. 
Every year, nearly 180,000 new prostate cancer cases are 
diagnosed, and prostate cancers in U.S annually cause about 
37,000 deaths. The developed cancers may spread to the 
lymph nodes or bones causing persistent and increasing 
pain, abnormal function, and death. The detection and 
treatment of early small prostate cancers are most important 
to prevent death attributable to prostate cancer. Current 
noninvasive approaches to detect the early prostate include 
the ultrasound, MRI and CT imaging, which have poor 
spatial resolution and contrast. Other means, such as needle 
biopsy, are invasive. Optical image as disclosed here can be 
used to image tumors in prostate. 

Optical tomography is being developed as a noninvasive 
method that uses nonionizing near-infrared (NIR) light in 
the 700-1500 nm range to obtain images of the interior of 
the breast and prostate. Tissues scatter light strongly, so a 
direct shadow image of any tumor is generally blurred by 
scattered light. A technique, known as inverse image recon- 
struction (HR), may help circumvent the problem of scat- 
tering. An HR approach uses the knowledge of the charac- 
teristics of input light, measured distribution of light 
intensity that emerges from the illuminated breast and pros- 
tate, and a theoretical model that describes how light propa- 
gates through breast and prostate to construct an image of 
the interior of breast and prostate. Back scattering and 
transmission geometries are used for breast. For prostate, 
backscattering geometry is more suitable. The image recon- 
struction methods can also be used to image objects in 
hostile environments of smoke, cloud, fog, ocean, sea, and 
to locate corrosion under paint. 

Although the problem has received much attention lately, 
that development of optical tomography has been slow. One 
of main difficulties is lack of an adequate algorithm for 
inversion image reconstruction, which is able to provide a 
3D image in reasonable computing times. Recent algorithms 
and methods have been developed to solve the inverse 
problem in order to produce images of inhomogeneous 
medium, including finite-element solutions of the diffusion 
equation and iterative reconstruction techniques, modeling 
fitting, least- square-based and wavelet based conjugate- 
gradient -decent methods. Examples of references which 
disclose this technique include: H. B. Jiang et al, “Fre- 
quency-domain optical image reconstruction in turbid 
media: an experimental study of single-target detectability,” 
Appl. Opt. Vol. 36 52-63 (1997); H. B. O’Leary et al, 
“Experimental image of heterogeneous turbid media by 
frequency -domain diffusing-photon tomography,” Opt. Lett. 
Vol. 20, 426-428 (1995); S. Fantini et al, “Assessment of the 
size, position, and optical properties of breast tumors in vivo 
by noninvasive optical methods,” Appl. Opt. Vol 36, 
170-179 (1997); W. Zhu et al “Iterative total least-squares 
image reconstruction algorithm for optical tomography by 
the conjugate gradient method,” J. Opt. Soc. Am. A Vol. 14 
799-807 (1997), all of which are incorporated herein by 
reference. These methods take significantly long computing 
time for obtaining a 3D image to be of use in clinical 
applications and other detections. These methods require a 
linearly or non-linearly inverse of a group of equations, 
which have huge unknown augments that equal to the 
number of voxels in a 3D volume (a voxel is a small volume 
unit in 3D volume, and corresponds to a pixel in 2D plane). 

Using a Fourier transform inverse procedure can greatly 
reduce computing time. Examples of references that disclose 
this technique include: X. D. Li et al, “Diffraction tomog- 
raphy for biomedical imaging with diffiise-photon density 
waves,” Opt. Lett. Vol. 22, 573-575 (1997); C. L. Matson et 
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al, “Analysis of the forward problem with difluse photon 
density waves in turbid media by use of a diffraction 
tomography model,” J. Opt. Soc. Am. A Vol. 16, 455-466 
(1999); C. L. Matson et al, “Backpropagation in turbid 
media”, J. Opt. Soc. Am. A Vol. 16, 1254-1265 (1999), all 5 
of which are incorporated herein by reference. In the Fourier 
transform procedure, the experimental setup should satisfy 
the requirement of spatial translation invariance, which 
restricts, up to now, use of a single laser source (a point 
source or a uniformly distributed plane source) with a 2D to 
plane of detectors in parallel (transmission or reflection) 
geometry. This type of experimental setup can acquire only 
a set of 2D data for continuous wave (CW) or frequency- 
domain tomography, which is generally not enough for 
reconstruction of a 3D image, resulting in uncertainty in the 15 
depth of the objects in 3D image. 

To overcome this difficulty of lack of enough data for 3D 
image in tomography using a Fourier procedure, we have in 
the past developed algorithms to acquire time-resolved 
optical signals, which provides an additional ID (at different 20 
times) of acquired data, so 3D image reconstruction can be 
performed. Examples of references which disclose this tech- 
nique include: R. R. Alfano et al: “Time-resolved diffusion 
tomographic 2D and 3D imaging in highly scattering turbid 
media,” U.S. Pat. No. 5,931,789, issued Aug. 3, 1999; U.S. 25 
Pat. No. 6,108,576, issued Aug. 22, 2000; W. Cai et al: 
“Optical tomographic image reconstruction from ultrafast 
time-sliced transmission measurements,” Appl. Optics Vol. 

38, 4237-4246 (1999); M. Xu et al, “Time-resolved Fourier 
optical diffuse tomography”, JOSAA Vol. 18 1535-1542 30 
(2001). Schotland and Markel developed inverse inversion 
algorithms using diffusion tomography based on the ana- 
lytical form of the Green’s function of frequency -domain 
diffusive waves, and point-like absorbers and scatterers. 
Examples of references which disclose this technique 35 
include: V. A. Markel, J. C. Schotland, “Inverse problem in 
optical diffusion tomography. I Fourier-Laplace inversion 
formulas”, J. Opt. Soc. Am. A Vol. 18, 1336 (2001); V. A. 
Markel, J. C. Schotland, “Inverse scattering for the diflusion 
equation with general boundary conditions”, Phys. Rev. E 40 
Vol. 64, 035601 (2001), all of which are incorporated herein 
by reference. 

From the viewpoint of data acquisition in parallel geom- 
etry, however, it is desirable to use a 2D array of laser 
sources, which can be formed by scanning a laser source 45 
through a 2D plane, and a 2D plane of detectors, such as a 
CCD camera or a CMOS camera. Each illumination of laser 
source produces a set of 2D data on the received detectors. 

For CW or frequency-modulated laser source, this arrange- 
ment can produce a set of (2D, 2D)=4D data in a relatively 50 
short acquisition time, with enough accuracy and at reason- 
able cost. When time-resolved technique is applied using a 
pulse laser source, a set of 5D data can be acquired. In these 
cases the inverse problem of 3D imaging is over-determined, 
rather than under-determined for the case of using a single 55 
CW or frequency domain sources, and thus produces a much 
more accurate 3D image. 

The key point is how to develop an algorithm, which is 
scientifically proper, and runs fast enough to produce a 3D 
image, so it can be realized for practical clinical applications 60 
and other field applications. 

SUMMARY OF THE INVENTION 

One object of the present invention is to teach and provide 65 
an inverse algorithm for fast three dimensional tomographic 
image reconstruction using a hybrid-dual -Fourier inverse 


4 

method suitable for experimental arrangements of multiple 
sources and multiple detectors in parallel (transmission or 
backscattering) geometries for either CW, frequency-do- 
main, and/or time-resolved approaches. 

Another object of the present invention is to teach and 
provide an inverse algorithm for three dimensional tomo- 
graphic image reconstruction using a hybrid-dual -Fourier 
inverse method based on experimental arrangements of 
multiple sources and multiple detectors in cylindrical geom- 
etries for either CW, frequency-domain, and/or time-re- 
solved approaches. 

A further object of the present invention is to teach and 
provide a novel hybrid-dual-Fourier mathematical method 
as an extension of the standard Fourier deconvolution 
method, applied to any kind of N-dimensional dual decon- 
volution problems where measurement data depend on two 
variables, and weight function satisfy the condition of trans- 
lation invariance for each variable in a M-dimensional 
subspace. 

Yet another object of the present invention is to teach and 
provide an accurate analytical solution of the Boltzmann 
photon transport equation in a uniform medium to serve as 
background Green’s function in the forward physical model 
for the tomography method of the present invention. 

Still another object of the present invention is to teach and 
provide a tomographic method using laser sources with 
different wavelengths for producing an internal map of a 
specific material structure in a turbid medium. 

One other object of the present invention is to teach and 
provide experimental designs for using hybrid-dual-Fourier 
tomography for detecting cancer and to develop an optical 
tomography imaging system. 

Additional objects, as well as features and advantages, of 
the present invention will be set forth in part in the descrip- 
tion which follows, and in part will be obvious from 
description or may be learned by practice of the invention. 

In accordance with one embodiment of the present inven- 
tion, a method for imaging an object in a turbid medium 
comprises the steps of: 

(a) directing an incident light from a source onto the 
turbid medium to obtain a plurality of emergent waves 
from the turbid medium; 

(b) determining the intensity data of at least part of the 
emergent waves by a plurality of detectors; 

(c) repeating the steps of a) and b) by placing the source 
of incident light at different positions until data acqui- 
sition is substantially complete; and 

(d) processing the intensity data by using an image 
reconstruction algorithm including a forward physical 
model and an inverse algorithm, to inversely construct 
a three dimensional image of the object in the turbid 
medium, where the inverse algorithm is a hybrid dual 
Fourier tomographic algorithm. 

In accordance with another embodiment of the present 
invention, a system for imaging an object in a turbid medium 
comprises: 

(a) a source for directing an incident wave onto the turbid 
medium to obtain a plurality of emergent waves from 
the turbid medium; 

(b) a plurality of detectors disposed along the propagation 
paths of at least part of the emergent waves for deter- 
mining the intensity data of the emergent waves; and 

(c) a data processor connected to the detectors to process 
the obtained intensity data and produce a three dimen- 
sional image of the object in the turbid media, where 
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the data processor is programmed to execute an inverse 
algorithm based on a forward physical model to process 
the intensity. 

The present invention is directly related to an optical 
tomographic method for imaging hidden objects in highly 5 
scattering turbid media. In one aspect of the present inven- 
tion, a method for imaging objects in a highly scattering 
turbid medium in parallel geometry includes the steps of: 
using a light source in visible and/or infrared spectral region, 
step by step, scanning through a two dimensional (2D) array, 
to illuminate a highly scattering medium; in each scanning 
step, to acquire signals of transmitted or backscattered light 
emergent from the medium received by a two dimensional 
(2D) array of detectors, such as CCD camera or CMOS 
camera; and applying a novel hybrid-dual Fourier inverse 
algorithm to form a three dimensional image of the objects 
in the highly scattering turbid medium. FIG. 1 schematically 
shows the experimental setup in parallel geometry used for 
the teaching presented here. ^ 

Preferably, a novel hybrid-dual-Fourier inverse algorithm 
is developed based on the present inventor’s discovery, that 
by performing a dual 2D Fourier transform upon both 
arguments related to the source positions and the detector 
positions, and using a hybrid linear transform, 3D tomog- 2 5 
raphy can be realized for use of multiple sources and 
multiple detectors. The main concept is as follows. A linear 
forward model describing light migration through an inho- 
mogeneous scattering medium with parallel geometry, based 
on the Bom approximation, either for CW, frequency 30 
domain, and/or time-resolved approaches, can be written as: 


following, we make a dual 2D Fourier transform 

Jd r /l r ^e z qsrs e z q<fr<f on equation (1), and obtain that 

s,z*z s )=J dzW{ q d , q s ,z,z d ,z s )X{ q^q s ,z), (2) 

where Y, X, and W are the corresponding Fourier space 
quantities in equation (1). 

Equation (2) seems most difficult to be used for perform- 
ing the Fourier inverse reconstruction because the arguments 
of X are different from that of Y and W. To remove this 
complexity, we perform a linear hybrid transform of the 
detector’s and source’s spatial frequency coordinates: 

u = qj+<is 


v = c id 


that leads to the following formula: 

Y( u , v ,z d ,z s )=J dzW(u , v ,z,z d , z s )Jf( u ,z), 


(3) 


(4) 


Y( r d , r s ,z d ,z s )=J d r dzW( r d - r , r s - r ,z,z d ,z s ) X 
( r ,z\ 


( 1 ) 


35 


where R =( r ,z) denotes the position of a voxel inside 
turbid medium; r is (x, y) coordinates; R 5 =( r 5 ,z^)denotes 

the position of a source; R d ={ r d ,zJ) denotes the position of 

— > — > 40 _ _ 

a detector. In equation (1),Y( r d , r s ,z d ,z s ) is the measured ID unknown value of X( u ,z) from known 2D data of Y( u , 


where Y, X, and W are, respectively, Y, X, and W as 
functions of u and v . 

FIG. 2 schematically explains the linear hybrid transform 
in equation (3), using an example of 6x6 lattices, from (q^, 
q^) coordinates to (u, v) coordinates. Note that the periodic 
property of lattices in the Fourier space is used, for example, 
Y(u=2, v=4)=Y(q rf =3, q^=5) as shown in FIG. 2. This figure 
shows that Y and W at each node in (u, v) coordinates can 
be obtained, respectively, from Y and W at the correspond- 
ing node in (q rf , qj coordinates without any algebraic 
manipulation. 

The hybrid transform equation (3) is a key, which makes 
the inverse reconstruction much easier to be performed. For 

each value u , equation (4) leads to an over-determining ID 
problem for inverse reconstruction, namely, to determine a 


change of light intensity, which incident from a source at \i s 

and received by a detector at R d . The word “change” refers 
to the difference in intensity compared to that received by 45 
the same detector, from the same source, but light passing 

through a homogeneous background medium; X( r ,z) is the 
change of the optical parameters inside turbid medium, in 
particular, it means change of the absorption coefficient j x a 
and the reduced scattering coefficient \ x/ in diffusion tomog- 
raphy. W(r rf -r,r 5 -r ,z,z d ,z s ) is the weight function, 

which is function of r r and r s - r at (x, y) plane, 
because of parallel geometry, and the translation invariance 
of the Green’s function in a homogeneous background 
medium. Here, we do not specify what form of the weight 
function; it can be an expression of the diffusion forward 
model, either for CW, Frequency domain, and/or time- 
resolved cases, or one based on the cumulant analytical 
solution of the radiative transfer equation we recently devel- 
oped. 

The inverse problem is to determine value of X from 
known measured data Y. The common understanding is that 
since the weight function now is related to three positions: 

r d , r s , and r , translation invariance cannot be simulta- 
neously satisfied, hence, it is difficult to performing Fourier 


v ) for each u . Detail of this ID procedure is described in 
the section of “detailed description of preferred embodi- 
ments.” This task is much easier than direct inversion of 

equation (1), which is a 3D inverse problem. After X( u ,z) 

for all u are obtained, a 2D inverse Fourier transform 

produces X( r ,z), which is the 3D image of optical param- 
50 eters of the body. 

Using this algorithm to perform an inverse reconstruction 
of 3D image of breast tissue with enough fine resolution (for 
example, 32x32x20 voxels) only takes a few minutes on a 
personal computer. 

55 Preferably, the hybrid-dual-Fourier inversion method can 
be used for turbid medium of substantially the cylindrical 
geometry, with an arbitrary shape of the (x, y) cross section, 
for 3D tomographic image reconstruction. FIG. 3 schemati- 
cally shows the experimental setup in the cylindrical geom- 
60 etry. Under this geometry, an algorithm using single-Fourier 
inversion was developed. Examples of references which 
disclose this technique include: Cai et al. “Three dimen- 
sional image reconstruction in high scattering turbid media,” 
SPIE 2979, p24 1-248 (1997), which are incorporated herein 
65 by reference. This algorithm limits to apply to the cases that 
the sources and the detectors are located at the same z plane, 
which restricts acquiring data. The following invention 
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provides a hybrid-dual -Fourier inverse approach for cylin- 
der geometry to remove the above-mentioned limitation, so 
more data can be acquired for 3D tomography. The linear 
forward model in the cylinder geometry is given by 

Y( r d , r s ,z d ,z s )=\d r dzW(r d , r s , r ;z d -z,z s -z)X 

Gix ( 5 ) 

where W( r d , r r ;z rf -z,z 5 -z) is the weight function, 
which is function of z^-z and z^-z because of cylinder 10 
geometry (assuming infinite z length) and the translation 
invariance of the Green’s function in a homogeneous back- 
ground medium. 

We make a dual ID (along z direction) Fourier transform 
fdzjlz s e iqjZd e iq ^ s on equation (5), and obtain that 15 

Yiq d ,q s , r d ,r s )=\ dzW(q d ,q s , r , r d , r s )X{q^q s , r ), (6) 

where Y, X, and W are the corresponding Fourier space 
quantities in equation (5). We further perform a linear hybrid 20 
transform (ID) of the detector’s and source’s spatial fre- 
quency coordinates: 

u=q d +q s 

v=qj-q„ (7) 25 


that leads to: 


Y{u,v, r d , r 5 )=J d r W(u, v, r d , r s ; r )X(u, r ), (8) 

where Y, X, and W are, respectively, Y, X, and W as 30 
functions of u and v. For each value of u, the above- 
mentioned equation leads to a over-determining 2D problem 
for inverse reconstruction, namely, to determine a 2D 

unknown value of X(u, r ) from known 3D data of Y(u,v, 35 


r rf ,rj for each u. This 3D-2D determination enhances 
accuracy of 3D image comparing to 2D — 2D determination 

in the single-Fourier transform inversion. After X(u, r ) for 
all u are obtained, a ID inverse Fourier transform produces 40 


X( r ,z), which is the 3D image of optical parameters. 

Preferably, the formula of the hybrid-dual -Fourier 
approach, equation (1) through equation (4), can be regarded 
a pure mathematical method to solve a N-dimensional dual 45 
deconvolution problem as an extension of the standard 
Fourier deconvolution method. In a N dimensional space, a 
deconvolution problem is defined by equation (1), where 

R =( r ,z) with r in a M subspace and z in a N-M subspace; 

Y( r d , r s ,z d ,z s ) is the measurement data which depend on 

both R^=( r s ,z s ) and R rf =( r d ,z d ). X( r ,z) is the quantity 

that should be determined by deconvolution. W( r r , 

r s - r ,z,z d ,z s ) is the weight function, which is function of 


50 


55 


r r and r s - r at M-dimensional sub space. 

The approach described from equation (1) through equa- 
tion (4), hence, can be applied to many other applications 
that lead to the form of equation (1 ). In this form, the sources 60 
can be light, X-ray, microwave, sound, electrons, particles, 
mechanical vibration, etc; the detectors can be any type of 
sensors for receiving light, X-ray, microwave, sound, elec- 
tricity, mechanical signals, etc; the “space” can be positions, 
times, wavelength spectrum, vibration modes, etc. Several 65 
examples are as follows: (1) using electrons from a linear 
electron accelerator through a cargo to detect merchandise 


inside cargo in the custom; (2) using sound to detect mine 
vertical distribution under ground; (3) using pressure vibra- 
tion on the surface of a material to detect the elastic 
coefficients inside the body. In the abovementioned and 
other examples, the source (electrons, sound, vibration, etc.) 
can be scanned on a 2D surface, and sensors can be arranged 
on the transmission and/or back-reflect 2D surface. The 
novel hybrid-dual -Fourier inverse algorithm can be applied 
in these cases for fast 3D imaging. 

Preferably, an accurate analytical solution of the Boltz- 
mann photon transport equation in an infinite uniform 
medium, first derived by the inventors, is combined to the 
above-mentioned inverse algorithm to provide more accu- 
rate forward model than the diffusion forward model. 
Examples of references which disclose this technique 
include: R. R. Alfano et al., “Time-resolved optical back- 
scattering tomographic image reconstruction in scattering 
media”, U.S. Pat. No. 6,205,353 issued Mar. 20, 2001; W. 
Cai, et al., “Cumulant solution of the elastic Boltzmann 
transport equation in an infinite uniform medium”, Phys. 
Rev. E. 61 3871 (2000); W. Cai, et al., “Analytical solution 
of the elastic Boltzmann transport Equation in an infinite 
uniform medium using cumulant expansion”, J. Phys. 
Chem. B104 3996 (2000); W. Cai, et al., “Analytical solu- 
tion of the polarized photon transport equation in an infinite 
uniform medium using cumulant expansion,” Phys. Rev. E 
63 016606 (2001); M. Xu et al., “Photon-transport forward 
model for imaging in turbid media,” Opt. Lett. 26, 
1066-1068 (2001), which are incorporated herein by refer- 
ence. The detailed description is given in the section of 
“detailed description of preferred embodiments.” 

In addition, this method can be used to determine the local 
material structure by distinguishing different values of opti- 
cal parameters obtained by using different light wave- 
lengths. Water, blood, and fat have different scattering and 
absorption parameters at different wavelengths in NIR 
region. Other biological materials, such as cancer, precan- 
cerous, and benign tissue will have different value of these 
optical parameters (scattering and absorption). For example, 
assume that cancer has the absorption and scattering param- 
eters at a wavelength X x different from that at a wavelength 
X 0 . When two sources are used having respective wave- 
lengths X 0 and X 1? where X 0 is a non-characteristic wave- 
length, the difference of their absorption coefficients p^(r, 
Xi)-|x a (r,Xo) and the scattering coefficients p^(r,Xi)-p^(r,X 0 ) 
can be obtained by inverse computation. This process pro- 
vides a significantly clearer image map of fat location by 
eliminating the background values. This procedure can yield 
maps of water, fat, blood, and calcification, even possibly 
cancer, using different X. 

Other objects and features of the present invention will 
become apparent from the following detailed description 
considered in conjunction with the accompanying drawings. 
It is to be understood, however, that the drawings are 
designed solely for purposes of illustration and not as a 
definition of the limits of the invention, for which reference 
should be made to the appended claims. It should be further 
understood that the drawings are not necessarily drawn to 
scale and that, unless otherwise indicated, they are merely 
intended to conceptually illustrate the structures and proce- 
dures described herein. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

In the drawings: 

FIG. 1 is a simplified schematic view of device for 
detecting breast cancer in parallel geometry using the 5 
hybrid-dual -Fourier tomographic algorithm of the present 
invention. 

FIG. 2 is a diagram for explaining the linear hybrid 
transform using an example of 6x6 lattices from (q^, qj 
coordinates to (u, v) coordinates of the present invention, to 

FIG. 3 is a simplified schematic view of device for 
detecting breast cancer in cylinder geometry using the 
hybrid-dual -Fourier tomographic algorithm of the present 
invention. 

FIG. 4 is a block/flow diagram of an optical tomography 15 
system/process in accordance with an embodiment of the 
present invention. 

FIGS. 5a, 5b, and 5c are simplified schematic views of 
devices for detecting breast cancer using the hybrid-dual - 
Fourier tomography method in parallel geometry of the 20 
present invention. 

FIG. 6 is a diagram of transillumination images of a 5 mm 
thick human breast tissue sample comprising adipose and 
fibrous regions obtained using light of different wave- 
lengths: (a) 1225 nm, (b) 1235 nm, (c) 1255 nm, and (d) 25 
1300 nm from a Cr:forsterite laser. This figure is taken from 
the reference: S. K. Gayen et al., “Near-infrared laser 
spectroscopic imaging: a step towards diagnostic optical 
imaging of human tissues”, Lasers in the Life Sciences Vol. 

8 187 (1999). 30 

FIG. 7 is a schematic diagram illustrating image maps of 
key components of breast tissue using different wavelengths. 

FIG. 8 is a diagram of comparative 3D images of a hidden 
absorbing object located at position (15, 15, 10) inside of a 
turbid medium divided into 32x32x20 voxels reconstructed 35 
by hybrid dual Fourier tomography in transmission parallel 
geometry based on the diffusion forward model. 

DETAILED DESCRIPTION OF THE 
PRESENTLY PREFERRED EMBODIMENTS 4 o 

The present invention is directed to novel optical tomo- 
graphic system and method for imaging hidden objects in 
highly scattering turbid media. Referring now to FIG. 4, a 
block diagram illustrates an optical tomography system in 45 
accordance with one aspect of the present invention. It is to 
be understood that the block diagram depicted in FIG. 4 may 
also be considered as a flow diagram of a method for 
imaging objects in turbid media in accordance with the 
present invention. The system 10 includes an illumination 50 
source 12, step by step, shining (direct scanning or using 
optical fiber) through a two-dimensional array on the plane 
surface of the turbid medium 30 in parallel geometry, or 
through around the cylinder surface of the turbid medium 30 
in cylinder geometry, and illuminating the turbid medium 55 
30. The illumination source 12 is a laser that emits a 
continuous wave, or a frequency-modulated light, or 
ultrashort light pulses (e.g., fsec, psec, and nsec pulses) 
having wavelengths in the range of about 700 to 1 500 nm so 
as to obtain deep penetration of the turbid medium 30 (such 60 
as breast, prostate, brain tissue, and cloud etc.). The laser 
source may include any conventional laser such as a semi- 
conductor laser, a Ti: Sapphire laser, a Cr 4+ Forsterite laser, 
a Cr 4+ YAG lasers, and a Cr 4+ —Ca 2 Ge0 3 (CUNYITE), a 
Nd:YAG laser. 65 

A plurality of detectors 14 located at the transmitted plane 
surface or the backscattering plane surface of the turbid 


10 

medium in parallel geometry, or located around cylinder 
surface of the turbid medium in cylinder geometry, are 
provided for acquiring signals of scattered light emergent 
from the turbid medium 30 for each shine of laser source. 
The detectors 14 are implemented CCD (charge coupled 
device) system or a group of fiber-detectors, and in the case 
of time-resolved measurement a time gating Kerr or inten- 
sified CCD is used for detecting pico-second time slicing 
signals. The light signals (“intensity data”), which are 
received by detectors 14, are intensity as functions of the 
position of the source 12 and detector 14, as well as the 
injecting direction of the source 12 and the receiving direc- 
tion of the detector 14. 

The intensity data which are detected and collected are 
processed via an inverse computation module 16 using a 
novel fast hybrid dual Fourier reconstructing algorithm to 
produce a three-dimensional image map of the internal 
structure of the turbid medium 30. The reconstruction algo- 
rithm (which is utilized by the inverse computation module 
16) includes a forward physical model 20. The forward 
model 20 (which is discussed in further detail below) 
describes photon migration (light propagation) in the turbid 
medium in accordance with optical parameters characteristic 
of a turbid medium: scattering rate, absorption rate, and, 
possible, differential angular scattering rate. The forward 
model 20 is based on an analytical solution 22 to the 
Boltzmann photon transport equation, or its diffusion 
approximation. Specifically, the analytical solution 22 com- 
prises a cumulate solution of the Boltzmann photon trans- 
port equation or a diflusive solution in an infinite uniform 
medium and a corresponding solution in a slab uniform 
medium, by adding virtual sources. The analytical solution 
22 serves as the background Green’s function of the forward 
physical model 20 for the present tomographic method. 

An inverse algorithm module 18, which employs a novel 
hybrid-dual -Fourier inverse algorithm, unique to the present 
imaging method, generates an internal map of the turbid 
medium by reconstructing the turbid medium structure. The 
inverse process is discussed below in further detail. 

The reconstruction algorithm of the present invention 
includes a regularization module 24 that provides suitable 
regularization parameters for use by the inverse algorithm 
module 18. Conventional methods such as the L-curve 
method disclosed in “The Truncated SVD as a Method of 
Regularization,” by Hansen, BIT, 17, 354-553, 1987, and 
the generalized cross validation (GCV) method disclosed in 
“Generalized Cross-Validation as a Method for Choosing a 
Good Ridge Parameter”, by Golub et al., Technometrics, 21, 
p. 215-223 (1979), may be used in the regularization 
module 24 for providing suitable regularization parameters. 
These methods disclosed in these references are incorpo- 
rated herein by reference. 

The system 10 may also include a knowledge catalog 
system 26 for building a relationship between different 
tissue structures and their corresponding optical parameters 
at different wavelengths of light source. The catalog system 
26 is utilized by the inverse computation module 16 to 
determine the local tissue structure and refine the corre- 
sponding optical parameters at a position. This system 26 
can be utilized to determine the local material structure by 
distinguishing or determining the local material structure 
from the local optical parameters. 

The reconstruction algorithm of the system 10 also 
includes an image graphic display module 28 for generating 
and displaying 3-D reconstructed images. 
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Forward Physical Model 
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It is to be understood that the present system and method 
is preferably implemented on a fast speed computer, for 
example, PC or Silicon Graphic (SGI), for fast numerical 
computation and graphic display. 

It is to be further understood that the present system and 5 
method may be used to image various highly scattering 
turbid media such as biological plant tissue, animal tissue, 
and human tissue. With regard to human tissue, for example, 
the present invention can be utilized to image breasts, brain, 
prostate, arteries, liver, kidney, bones in joints, calcification 
regions, and arthritis, fingers, arms, legs, feet, etc. The turbid 
media, inside or through which the objects may be imaged, 
also includes cloud, fog, smog, dust, smoke, etc, as well as 
defects in semiconductors, ceramics, dielectrics, corrosion 15 
under paint. 

Inverse Algorithm 

A novel hybrid dual Fourier inverse reconstruction algo- 20 
rithm is developed based on the present inventor’s discov- 
ery, that by performing a dual Fourier transform upon both 
arguments related to the source positions and the detector 
positions, and using a linear hybrid transform, 3D tomog- 25 
raphy can be realized for case of multiple sources and 
multiple detectors. 

The formula for the parallel geometry is presented in 
equation (1) to equation (4) in the section “summary of the 
invention.” FIG. 2 schematically explains the linear hybrid 30 
transform in equation (3), using an example of 6x6 lattices 
from (q rf , qj coordinates to (u, v) coordinates. Note that the 
periodic property of lattices in the Fourier space is used, for 
example, Y(u=2, v=4)=Y(q i =3, q/=5)as shown in FIG. 2. 
This figure shows that Y and W at each node in (u, v) 35 
coordinates can be obtained, respectively, from Y and W at 
the corresponding node in (q rf , qj coordinates without any 
algebraic manipulation. 

The inverse formula can be written as 

40 

X=Y t W[W t W+ AJ-\ For each a, (9) 

which inverses a matrix W r W with N z rank, where N z is 
number of 1 D division along z layer. In equation (9), A is the ^ 
regularization matrix. Examples of references which dis- 
close the regularization technique include: R. R. Alfano et al. 
“Time-resolved diffusion tomographic 2D and 3D imaging 
in highly scattering turbid media,” U.S. Pat. No. 5,931,789, 
issued Aug. 3, 1999; U.S. Pat. No. 6,108,576, issued Aug. 

22, 2000, all of which are incorporated herein by reference. 

After X( u ,z) for all u are obtained, a 2D inverse Fourier 

transform produces X( r ,z), which is the 3D image of 
optical parameters. By use of the abovementioned procedure 55 
of tomography, 3D imaging is realized for multiple-sources 
and multiple-detectors in the parallel geometry. Using this 
algorithm to perform an inverse reconstruction of 3D image 
of breast tissue with enough fine resolution (for example, 
32x32x20 voxels) only takes a few minutes on a personal 60 
computer. 

The formula for the cylinder geometry is presented in 
equation (5) to equation (8) in the section “summary of the 
invention”. In this case, a ID dual Fourier transform is 
performed along z direction, and the matrix W r W in equa- 65 
tion (9), which should be inversed, has N xy rank, where N xy 
is number of 2D division at x-y plane. 


Forward Model Based on the Solution of the Radiative 
Transfer Equation 

The following discussion provides the theoretical basis 
for the forward model of the present invention. The structure 
of a highly scattering turbid medium can be characterized by 
the following optical parameters: |i 5 (r) the scattering rate; 
p, a (r) the absorption rate; and |i 5 (r)P(s',s,r) the differential 
angular scattering rate. Hereafter r denotes a 3D vector. 
These parameters are position dependent, and represent the 
non-uniform structure of the highly scattering turbid 
medium. The values of these optical parameters vary using 
light sources with different wavelengths, X. For instance, the 
absorption rate, p^(r) will vary with the wavelength because 
the absorption peak appears when the wavelength matches 
the difference of the energy levels of a specific molecular 
structure. In addition, the scattering rate, fi/r), and the 
differential angular scattering rate, p^(r)P(s , ,s,r) vary with the 
wavelength because these rates are related to R/X, where R 
is the average radius of the scatterer. 

The photon propagation in a medium is described by the 
photon distribution function, I(r,s,t), namely, the photon 
intensity in a unit of solid angle as functions of time t, 
position r, and direction s. The mathematical equation gov- 
erning photon propagation is the well-known Boltzmann 
radiative transfer equation: 

dl(r,s, t)/dt+cs •'V r I(r,s, t)+\i a {r)I{r,s, ?)=}i s (r) J P(s,s\ r ) [/ 

{r,s't)-I{r,sj)]ds'+b(r-r<y)b{s-s^) 6(?-0) (10) 

It is difficult to directly solve the above radiative transfer 
equation. A perturbation method is used which designates 
the photon distribution function in a uniform background 
medium as the zero -order approximation. This method des- 
ignates, as the first-order perturbation, the change of the 
photon distribution function due to the change of optical 
parameters compared to that in the uniform background 
medium. The change of scattering and absorption param- 
eters are defined as follows: 

A ii s (r)=ii s (r)-ii s ^; 

and 

AQi^P] ( s \ s,r)=[L s (r)P(s\s,r)-\i s ( - 0 ^P^ 0 \s ’,s); (1 1) 

where the quantities with super index (0) are the optical 
parameters in a uniform background medium (a medium 
without hidden objects). By expanding A[|x 5 .P](s',s,r) in Leg- 
endre polynomials, we get: 


A[ji s P]{s'j,r) = ^-Y J ^t i s(r)^a l (r)P l [cos(s , s)] 


with normalization of Aa 0 (r)=l. The corresponding Leg- 
endre coefficients, Aji/rJAa^r) can also serve as optical 
parameters. The following equation based on the standard 
Bom approximation method represents our forward model: 

M{r d ,s d , t\r s ,s s )=\ dt’\ dr\ ds 'f°\r d , s d , t-tV, s'){/A [M 3 ] 

( s',s , r)fi°\r,s, r1r s ,5 5 )d , 5-/A i u s (r)+A|i a (r)]/ 0) (?;5 ' 

tK*s)} ( 13 ) 

where AI (r rf ,s rf ,tlr^,sj is the change in light intensity 
received by a detector located at r d , along the direction s d , 
and at time t, which is injected from a source located at r s , 
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along a direction of s s , at time t=0. The word “change” refers 
to the difference in intensity compared to that received by 
the same detector, from the same source, but light passing 
through a uniform background medium (i.e., a medium 
without hidden objects). The term I (0) (r 2 ,s 2 ,tlr 1 ,s 1 ) is the 5 
intensity of light located at r 2 along the direction s 2 an at time 
t, which is injected from a position r 1 along a direction of s x 
at time t =0 migrating in a uniform background medium. 
Examples of references which disclose technique for obtain- 
ing I (0) (r^s^tlr^sj in an infinite uniform medium include: to 
“W. Cai, et al., “Cumulant solution of the elastic Boltzmann 
transport equation in an infinite uniform medium”, Phys. 
Rev. E. 61 3871 (2000), W. Cai, et al., “Analytical solution 
of the elastic Boltzmann transport Equation in an infinite 
uniform medium using cumulant expansion”, J. Phys. 15 
Chem. B104 3996 (2000), W. Cai, et al., “Analytical solu- 
tion of the polarized photon transport equation in an infinite 
uniform medium using cumulant expansion,” Phys. Rev. E 
63 016606 (2001), which are incorporated herein by refer- 
ence. Examples of references which disclose technique for 20 
obtaining (r 2 ,s 2 ,tlr 1 ,s 1 ) in an semi-infinite and slab 
shaped uniform medium include: R. R. Alfano et al., “Time- 
resolved optical backscattering tomographic image recon- 
struction in scattering media”, U.S. Pat. No. 6,205,353 
issued Mar. 20, 2001, which is incorporated herein by 25 
reference. 

We expand the background Green’s function in spherical 
harmonics: 


X YMYLV) = P/[cos(5-i')], 

the analytical integration over s and s’ in Eq. (13) can be 
performed. For time resolved data, the contribution from an 
absorbing object located at r k is given by 


A I{r d , s d , r s , s s , t\r k ) = -Afi a (r k )6V k 


dt ' 


Jo 


(15) 


L 

-Am (A? ^s? •$$? ? fy, S d , t O 

1=0 m 


where bV k is the volume of k^ voxel, and L is the cut-off 
value in the Legendre expansion in Eq. (15). The contribu- 
tion from a scattering object located at r k is given by 


A I(r d ,s d , r s ,s s , l\r k ) = 



4 n 

(2/TTj 


A^ s (r*) 1 - 


21 + 1 


( 0 ) Aa/(^) 
21 + 1 


(16) 


2 A/nfa? r s , t’)C* lm {r k , r d , s d , t-t’) 


I {0) (r, s, t'\r s , s s ) = ^jA lm (r, r s , s s , + B lm (r, r s , s s , Y)Y { ^\s), 

l,m 

I {0) (r d , s d , t- Y \r, s') = ^ Cimir, r d , s d , t-Y)Y ( ^{s , ) + Di m {r, r d , s d , t - Y )Y$ (s' ), 


where 


Y$(G, <p) = J D i" ,) (cos0)cos (m<p) 


For Frequency domain (or CW) data, the contribution 
from an absorbing object located at r k is given by 

40 


and 


A I(r d ,s d , r s ,s s , a)\r k ) = 


(17) 


Yj°\0, <p) = (cos0)sin(m0), 


with 


45 

L 

-A fi a (r k )6V k ^^ s *' ^ c lm( r k, r d , s d , w), 

i=o m 


P ( / " ,) (cos G) 


50 


the associated Legendre function. The spherical transform is 
performed using a fast Fourier transform for the integral 
over (|), and a Clenshaw-Curtis quadrature for the integral 55 
over 0 . 

Using the orthogonality relation of the spherical function and the contribution from a scattering object located at r k is 
and the addition theorem: given by 


M{r d , s d , r s , io\ r k ) = - SV k 


L 


4 n 


/ i2l+\ 
1 = 1 


A^(r*) 1 


a) 

' 21 + 


-Bs 


(0) Aa f (^) l 


2/ + 1 


( 18 ) 


^ A fair k , r s ,s s , a>)c] m {r k , r d ,s d , w). 
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By replacing Y=[AI/I (0) ] by -ln[I/I (0) ]), our model, to 
some extent, automatically includes higher order non-linear 
contribution. This procedure is usually called the Rytov 
approximation. 

A simplified photon-transport forward model has been 5 
developed. An example of references which disclose tech- 
nique to obtaining this forward model includes: M. Xu et al: 
“Photon-transport forward model for imaging in turbid 
media,” Opt. Lett. 26, 1066-1068 (2001), which is incor- 
porated herein by reference. 10 

Forward Model Based on the Solution of the Diflusion 
Equation 

By making Legendre expansion of the Boltzmann equa- 
tion (10) and cut at the lowest order, a diffusion equation is 15 
obtained as an approximate equation of radiative transfer: 


+ cn„(r) - V (cDMV)W, r) = S(r, t ) 


Where N(r,t) is the photon density, D(r) is the diffusive 
coefficient, and S(r,t) is the source distribution. The corre- 
sponding equation for cases of steady state and frequency- 25 
domain can easily derived from equation (19). 

Under first-order perturbation, the change of the photon 
density is determined from the change of optical parameters 
compared to that in the uniform background medium. The 
change of scattering and absorption parameters are defined 30 
as follows: 

AD(r)=D(r)-D (0) ; 

Afi a (r)=ii a (r)-iiJ 0) . (20) 35 

The corresponding forward model for the time-resolved 
case is given by 


A N(r d , r s , t) = J drcAju a (r)J drG°(r d , r, t - r)G°(r s , r, r) - 
J drAD(r)j^ df V G°(r d , r, t - r) • V G°(r s , r, r) 


( 21 ) 


40 


45 


The time-resolved Green’s function in an infinite uniform 
background medium is: 


G°(r, t) = 


(4^D ( °)cf) 3/2 eXP | 


(x 2 +y 2 +z 2 ) 
4 D^ct 


exp(-ju< 0 >r) 


( 22 ) 


50 


The corresponding forward model for the frequency- 
domain case is given by 

AN(r d ,r s ,(ti)=\ drcAju a (r)G°(r d , r,(ti)G°(r s , r,(o)-J dr AD 

(r)VG°(o,r,(o)-VG°(r 5 ,r,(o) (23) 

The frequency-domain Green’s function in an infinite 
uniform background medium is: 

G°(r,(o)=exp{ 1/2 Irl}/ (4 jiD (:o Vl) (24) 

Equations (22) and (24) can be extended to semi-infinite and 
slab medium by introducing the corresponding “image” 65 
sources. Equation (21) and (23) provide the needed weight 
function in equation (1) in the diffusion forward model. 
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Experimental Design for Image of Breast 


Referring now to FIG. 5, experimental devices are shown 
which may be utilized for detecting breast cancer using 
optical tomography method in parallel geometry of the 
present invention. As shown in FIG. 5(a), a sources-detector 
head 300, which includes a 2D array of sources and detec- 
tors, is fixed on a transparent plate 301. A medical doctor 
using a hand or other method (for example, moving the 
patient’s bed) can press the plate 301 against a patient’s 
breast to push the breast against the chest wall. Thereafter, 
a laser can be applied, and the detectors can then receive 
backscattered light signals. From these signals, through 
numerical computation by computer using the hybrid dual 
Fourier tomography algorithm of the present invention, a 
three-dimensional image of the entire breast can be recon- 
structed. Since breast is soft and flexible, it is possible to 
squeeze a breast to about 2 cm to 4 cm above the chest. In 
another embodiment as shown in FIG. 5(b), the breast may 
be squeezed between two parallel transparent plates 302. In 
addition, the embodiment shown in FIG. 5(c) can be used to 
detect a local breast region near the sources -detectors. The 
embodiment of FIG. 5u is similar to that shown in FIG. 5(a), 
but the source-detector head 300 and the plate 301 are 
smaller. By pushing successively upon different areas of the 
breast, a test of the entire breast can be completed. In order 
to reduce the clinic time, data acquisition can be performed 
during the visit with a patient. The image reconstruction then 
can be computed in parallel during the patient’s waiting 
time. If a near real-time image reconstruction can be real- 
ized, the doctor may also see an image of the local region of 
patient’s breast immediately after the laser beam is applied. 
Since only a local region of the breast is compressed (using 
the embodiment of FIG. 5c), it is possible to push down the 
local region of breast to 1 cm to 2 cm above the chest wall 
using a moderate pressure. The advantage for the embodi- 
ments of FIGS. 5(a) and 5(6) is that the image of whole 
breast can be reconstructed at one time, hence, the clinic is 
fast. The embodiment of FIG. 5(c), on the other hand, can 
enhance the image resolution and can reduce pain because 
only a local region of breast is tested at a particular time. 

Multiple Wavelengths 

Advantageously, this system can determine the local 
material structure by distinguishing different values of opti- 
cal parameters obtained using different light wavelengths. 
FIG. 6 shows the experimental images of transmission light 
from a human breast tissue sample, comprising adipose and 
fibrous regions, with light sources at different wavelengths. 
Publication which discloses this result includes S. K. Gayen 
et al., “Near-infrared laser spectroscopic imaging: a step 
towards diagnostic optical imaging of human tissues”, 
Lasers in the Life Sciences Vol. 8 187 (1999), which is 
incorporated herein by reference. The absorbing coefficient 
and the scattering coefficient through a type of tissue has 
different values with a laser source of different wavelengths. 
When two sources are used having wavelengths K 0 and X l5 
where k 0 is a non-characteristic wavelength, the difference 
of their absorption coefficients |i a (r,)q)-|i a (r,}q ) ) and the 
difference of scattering coefficients |i 5 (r,X 1 )-jq 5 (r,X 0 ) 
obtained by our inverse computation shows a more clear 
image map where the hidden object is located by eliminating 
the background values. This procedure of using different X 
can yield maps of water, blood, and calcification, cancer, 
precancerous and benign tissue. A schematic diagram for 
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using the different wavelengths in obtaining the internal 
maps of different components in a breast is shown in FIG. 

7. 


Test of Hybrid Dual Fourier Imaging Method 5 

Using Simulating Data 


As a proof of the concept of the hybrid-dual -Fourier 
tomographic algorithm, a 3D image reconstruction is per- 
formed from simulated data using the diffusion forward to 
model. 

3D Image in Parallel Geometry Using Hybrid-Dual -Fourier 
Tomographic Method 

A slab turbid medium with the optical parameters of 
breast, the transport mean free path 1,=1 mm, the absorption 
length 1^=300 mm, and thickness z rf =40 mm, is divided into 
20 layers. A CW laser shines, step by step, through a 32x32 
2D array on the z s =0 plane, each shining light passes through 
the medium and is received by a 32x32 2D array of detectors 2Q 
on the z^=40 mm plane. The medium is divided into 
32x32x20 voxels, each 3x3x2 mm 3 . A hidden object, with 
absorption length 1^=50 mm and volume 3x3x2 mm 3 , is 
located at position labeled (15, 15, 10). The tomographic 
images are shown in FIG. 8 using the new algorithm and 25 
hybrid transform. The number at FIG. Ha labels the z layers 
counting form source to the detector (each layer is separated 
by 2 nun). The FIG. Hb is the amplified figures of 9 th , 10^, 

1 1 th layers. The simulated results show that the obtained 
image has maximum value of absorption coefficient at the 3Q 
correct 3D position of the object (15,15,10). At the nearest 
neighbor in x-y lattice (15, 14, 10) or (15, 16, 10) etc., the 
values of absorption coefficient decrease about 20%, and 
then further decrease to about 50% at (15, 13, 10) etc., which 
indicates the resolution is about 6 mm in the transverse x-y 35 
plane. Comparing the values of the absorption coefficient at 
voxels (15, 15, 9) and (15,15,11) with (15, 15, 10), we see 
that they decrease about 20%, and further decrease to about 
50% at (15, 15, 8) etc. Since each layer is separated by 2 
mm, the resolution along z direction is about 6 mm. In 4Q 
general, the axial resolution (along z direction) is poorer 
than lateral resolution [at (x, y) plane]. In the parallel 
transmission geometry, two Green’s functions in the weight 
function compensate each other when the z position of the 
hidden object changes, which leads to a poor sensitivity of 45 
photon intensity to the z position of the object. 

The embodiments of the present invention described 
above are intended to be merely exemplary and those skilled 
in the art should be able to make numerous variations and 
modifications to it without departing from the spirit of the 5Q 
present invention. All such variations and modifications are 
intended to be within the scope of the present invention as 
defined in the appended claims. 


We claim: 55 

1. A method for producing a three-dimensional image of 
objects in a turbid medium from a hybrid dual Fourier 
transform, comprising: 

(a) measuring intensity data at multiple detectors from 

multiple sources in parallel geometry; 60 

(b) transforming data at positions of the detectors and the 
sources into Fourier coordinates (q^, qj to create a dual 
two-dimensional x-y spatial Fourier transform; 

(c) performing a linear hybrid transformation of the 
Fourier coordinates (q rf , qj of the multiple sources and 65 
multiple detectors in accordance with a predefined 
relationship; 
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(d) performing a one-dimensional inverse reconstruction 
for each first value of the predefined relationship; and 

(e) generating a fast two-dimensional inverse Fourier 
transform based on each first value to produce three- 
dimensional image of objects in the turbid medium. 

2. The method of claim 1, wherein said predefined rela- 
tionship comprises the following relationship: 


v = q<r<ls\ 

wherein u and v are coordinates, and q^, q^ are the Fourier 
coordinates of the multiple sources and multiple detectors, 
respectively. 

3. The method of claim 1, wherein said each first value is 
u of the following relationship: 


v = q dr q s, 

wherein u and v are coordinates, and q d and q^ are the 
Fourier coordinates of the multiple sources and mul- 
tiple detectors, respectively. 

4. The method of claim 1, wherein said predefined rela- 
tionship comprises the following relationship: 

u= q^q s 

v= qj-q s '> 

wherein u and v are coordinates, and q^ and q^ are the 
Fourier coordinates of the multiple sources and mul- 
tiple detectors, respectively. 

5. The method of claim 1, wherein said each first value is 
u of the following relationship: 

u= q^q s 


v=q d -q s ; 

wherein u and v are coordinates, and q^ and q^ are the 
Fourier coordinates of the multiple sources and mul- 
tiple detectors, respectively. 

6. The method of claim 1, wherein computational burdens 
of image reconstruction are reduced by combining a dual 
two-dimensional Fast Fourier Transform inversion with a 
one-dimensional matrix inversion. 

7. The method of claim 1, further comprising the step of: 

performing the linear hybrid transform to provide a fast 

three-dimensional Fourier image reconstruction for 
multiple sources and multiple detectors. 

8. The method of claim 1, wherein the multiple sources 
are selected from the group consisting of light, X-ray, 
micro-wave, sound, electrons, particles and mechanical 
vibration. 

9. The method of claim 1, wherein the detectors comprise 
sensors for detecting signals of one of light, sound, elec- 
tricity and mechanical waves. 

10. The method of claim 1, wherein the intensity data is 
a function of a coordinate space selected from the group 
consisting of position, time, wavelength spectrum and vibra- 
tion mode. 
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11. The method of claim 1, wherein a hybrid dual Fourier 
tomographic method is implemented for at least one of 
continuous- wave sources, frequency domain sources and 
ultra-fast pulse sources. 

12. The method of claim 1, wherein steps (b) thru (e) are 
solved based on N-dimensional dual de-convolution. 

13. The method of claim 1, wherein the turbid medium 
has a cylindrical surface, and the multiple sources and 
multiple detectors are arranged circumferentially with 
respect to the cylinderical surface. 

14. The method of claim 13, wherein the hybrid dual 
Fourier tomographic method for the cylindrical surface is 
performed in accordance with the following relationships 

Y( r d , r s ,z d ,z s )=\ d r dzW( r d , r s , r ;z d -z,z -z)X 
( r ,z); 

wherein W( r d , r r \z d -z,z s -z) is a weight function, 
which is function of z^-z and z s -z; 

Y(q d ,q s , r d ,r s )=\ dzW(q d ,q s , r , r d , r s )X(q d +q s , r ) 

wherein Y, X, and W are corresponding Fourier space 
quantities; 


u=q<&q s 

v=q d ~q s \ 

wherein u and v are coordinates, and q d , q^ are the Fourier 
5 coordinates of the source and detector, respectively; 
and 

Y{u,v, r d , r S )=J d r W(u, v, r d , r s ; r )X{u, r ); 

10 wherein Y, X, and W are, respectively, Y, X, and W as 
functions of u and v. 

15. The method of claim 1, wherein a forward physical 
model of the method is based on an analytical cumulant 
solution to a Boltzmann photon transport equation. 

15 16. The method of claim 15, wherein a background 

Green’s function of the method comprises the cumulant 
solution of the Boltzmann photon transport equation. 

17. The method of claim 15, wherein an accurate analyti- 
cal form of a solution of the Boltzmann equation is imple- 
20 mented to permit computation of the forward physical 
model. 



